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Abstract 

A matrix factorization problem is considered. The matrix to be 
factorized is algebraic, has dimension 2x2 and belongs to Moiseev's 
class. A new method of factorization is proposed. First, the matrix 
factorization problem is reduced to a Riemann-Hilbert problem us- 
ing the Hurd's method. Secondly, the Riemann-Hilbert problem is 
embedded into a family of Riemann-Hilbert problems indexed by a 
variable b taking values on a half-line. A linear ordinary differential 
equation (ODE1) with respect to b is derived. The coefficient of this 
equation remains unknown at this step. Finally, the coefficient of the 
ODE1 is computed. For this, it is proven that it obeys a non-linear 
ordinary differential equation (ODE2) on a half-line. Thus, the nu- 
merical procedure of matrix factorization becomes reduced to two runs 
of solving of ordinary differential equations on a half-line: first ODE2 
for the coefficient of ODE1, and then ODE1 for the unknown function. 
The efficiency of the new method is demonstrated on some examples. 



1 Introduction 



Many diffraction problems can be transformed into matrix factorization prob- 
lems pp. Typically, these diffraction problems are 2D problems with differ- 
ent boundaries occupying positive and negative parts of the x-axis. There 
emerges a known matrix G analytical in a thin strip going along the real axis 
of a complex variable k, and it is necessary to represent it as a product 

G{k) = U-\k)W{k), (1) 
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where W and U are matrices analytical and having no zeros of the deter- 
minant in the upper and lower half-plane, respectively. Also, both matrices 
should have algebraic growth in corresponding half-planes. 

In the scalar case, which can be considered as a degenerate case of 1 x 1 
matrix, the solution can be readily achieved by taking the logarithm of the 
matrix and performing the additive decomposition by means of Cauchy's 
integral and Sokhotsky's formula [lj. Returning to the matrices of order 
N > 1, this approach can be generalized for Moiseev's matrices [2] having 
form 

JV-l 

G=Y,9n(k)A n (k), (2) 

n=0 

where g n (k) are scalar function, and A is a polynomial matrix. In the simplest 
case of matrix A having distinct eigenvalues almost everywhere, A can be 
decomposed as 

A = TDT~ l , (3) 

where T is the matrix of the eigenvectors and D is a diagonal matrix com- 
posed of the eigenvalues. Both T and D are algebraic matrices, and one can 
introduce the Riemann surface TZ on which T and D are single-valued. Fur- 
ther, the matrix factorization problem becomes reduced to a scalar Riemann- 
Hilbert problem on TZ. This problem can be solved in terms of Abelian inte- 
grals with the help of Jacobi's inversion problem [3J. So, the solution of the 
problem of factorization of ([2]) is known at least formally, and it possibly can 
be used for practical needs. Some examples can be found e.g. in [3]. Simpler, 
but more popular cases [5j [6] can be described as particular cases of (J2]). 
Khrapkov's method [5] is rather simple and leads to straightforward compu- 
tations, but for a broad class of matrices it produces non-algebraic growth 
at infinity. Moiseev's method can be considered as a remedy enabling one to 
avoid this growth. Another technique to avoid the non-algebraic growth has 
been proposed in [TJ. This technique also includes some numerical stages. 
A review of the commutative factorization and a development of ideas of j6] 
can be found in [8]. 

If G cannot be represented as (|2J) then some numerical [9J or approximate 
(e.g. [10]) methods can be applied. 

In the current work we consider matrices of Moiseev's class (T2]) and de- 
velop a new technique which is arguably simpler in practical realization than 
the Moiseev-Zverovich or Daniele procedure. The new technique can be 
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applied only when the g n (k) in tfZ§ are algebraic functions. This is an impor- 
tant restriction, however in much of the practical situation this restriction 
is fulfilled. The new procedure comprises three steps. Firstly, the matrix 
factorization problem is reduced to a Riemann-Hilbert problem using the 
Hurd's method pTj . Namely, instead of studying the factors W and U we 
are studying only the factor U continued into the upper half-plane of k. The 
fc-plane is cut along half-lines connecting the branch points of G located in 
the upper half-plane, namely the points kj, with k = -Moo. A Hilbert prob- 
lem is formulated on the half-lines (kj, kj + zoo). As Hurd mentioned, the 
new problem can be simpler than the initial matrix factorization problem. 

At the second step the Riemann-Hilbert problem is embedded into a 
family of Riemann-Hilbert problems indexed by a variable b taking values 
on a half-line. Namely, for the whole family the coefficients Hj(k) remain the 
same, but the contours on which the functional equations should be fulfilled 
are changed from (kj, kj + ioo) to (kj + b, kj + ioo), where b is an imaginary 
number taking values from to -Moo. Thus, we can define the family of 
solutions U(b, k). The solution of the initial problem is denoted by U(0, k). 
A linear ordinary differential equation (ODE1) with respect to b is derived 
for U(b, k). The coefficient of this equation remains unknown on this step. 

Finally, step the coefficient of the ODE1 is computed. For this, it is proved 
that it obeys a non-linear ordinary differential equation (ODE2) on a half- 
line. Thus, the numerical procedure of matrix factorization becomes reduced 
to solving two ordinary differential equations on a half-line: first ODE2 for 
the coefficient of ODE1, and then ODE1 for the unknown function. 

Some numerical results are presented. Namely, we demonstrate that the 
new procedure applied to a matrix belonging to the Khrapkov's class is fac- 
torized exactly the same way as by the traditional Khrapkov's procedure. 
Moreover, we apply our method to the matrix emerging in [3]. 

2 Problem formulation and Hurd's procedure 

Let G(k) be an algebraic matrix N x N having no singularities and no zeros 
of determinant on the real axis and tending to the unit matrix / of dimension 
N x iV as | A; | — > oo. Our aim is to find the decomposition (JT|) valid in some 
strip |Im[A;]| < e with U having no singularities or zeros of the determinant 
in the lower half-plane and on the real axis, and W having no singularities 
or zeros of the determinant in the upper half-plane, maybe except several 
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points, where poles or zeros are allowed. We demand that the unknown 
functions U and W tend to I as \k\ — > oo. Some restrictions on G will be 
imposed below. 

Apply Hurd's procedure [TT] as follows. Let kj, j = 1, ... ,p be branch 
points of matrix G in the upper half-plane. Connect the points kj with zoo 
by the cuts 

Tj = (kj, kj + zoo) 

parallel to the imaginary axis. Let contours Tj do not pass through other 
branch points, poles or zeros of the determinant of G. Continue function 
U(k) into the upper half-plane cut along the lines Tj by the relation 

U{k) = W{k)G~ x {k). (4) 

Note that W is defined and regular in the upper half-plane and G is defined in 
the upper half-plane with the cuts Tj. Define by U(k + ) and U{k~) for k G Tj 
the values of U taken on the right and on the left shore of Tj, respectively 
(see Fig. [1]). Similarly, define the values G(k + ) and G(k~). For some k G Tj 

U(k + ) EE W(k)G'\k + ), 

U(k~) ee W(k)G-\k-) 
(note that W(k~) = W(k + ) = W(k)). Then, 

U{k + ) = U{k-)H J {k), (5) 

Hj(k) ee G(kr)G-\k + ), k G Tj (6) 

The set of equations (J6]) taken for j = l,...,p constitute the Riemann- 
Hilbert problem in Hurd'd formulation. It was Hurd's observation that this 
problem can be simpler than the initial matrix factorization problem. 

In our case each coefficient Hj(k) of the Riemann-Hilbert problem can 
be continued analytically into some half-strip Q + kj 

n = {Re[k] < e, lm[k] > 0}. 

Note that the point kj does not belong to Q + kj. Also Hj(k) are algebraic 
functions, so each of them can be continued onto some Riemann surface. 

Let all Hj(k) tend to I lm[k] — > oo (this restriction is fulfilled if G(k) — > I 
on all sheets). We forget about W(k) and look for U(k) on the complex 
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Figure 1: Cuts Vj 

plane k cut along the contours Tj having no singularities and no zeros of 
the determinant on the cut plane, obeying the problem ([5]). We recall that 
U(k) / as \k\ -> oo. 

The behavior of U(k) at the points kj will be specified below in such a 
way that the problem possesses a unique solution. 

We assume also that matrices Hj have distinct eigenvalues almost every- 
where. 

All restrictions described above correspond to a quite general matrix fac- 
torization problem and are easy to fulfil by, e.g. slight change of the contour 
position. If matrices Hj do not tend to / the method can be easily modified 
also. Here we are going to pose the most strong restriction: we assume that 
all branches of all matrices Hj(k) taken for arbitrary affix k commute with 
each other, i.e. for each k 

H n {k)H n {k)=H ]2 {k)H n {k), (7) 

where Hj 1 (k) and Hj 2 (k) are any possible continuations of Hj 1 and Hj 2 to k. 
The meaning of this restriction is discussed in the next section. 

3 Functional— commutative and branch— com- 
mutative matrices 

All existing analytical approaches to matrix factorization are available only 
for matrices admitting a commutative factorization, i.e. a representation of 
the form 

G{k) — U~ l {k)W(k) — W(k)U~ 1 (k). (8) 
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The theory of commutative matrix factorization starts from [12] where a con- 
cept of functional-commutative matrix has been introduced. A functional- 
commutative matrix is a matrix commuting with its singular integral. The 
property of functional-commutativity is not easy to check for an arbitrary 
matrix. That is why, in (13] we introduced branch-commutative matrices. 
Namely, an algebraic matrix G(k) is called branch-commutative if for any 
k the values Gj(k) corresponding to different branches of G commute with 
each other. 

To formulate the main result of [13] we need one more definition. A 
Riemann surface of an algebraic matrix is called balanced if each sheet of 
it can be reached from any fixed sheet only by bypassing the branch points 
located in the upper half-plane and only bypassing the branch points lying 
in the negative half-plane. Most of the known matrices arising in practical 
problems have balanced Riemann surfaces. 

The main result of [13J is as follows. If an algebraic matrix G with bal- 
anced Riemann surface admits commutative factorization then it is branch- 
commutative. Vice versa, a branch-commutative matrix can be represented 
in the form (J2J, and thus the Moiseev's method can be applied to it. 

Note that if matrix G is branch-commutative then the property of (J7J) 
for the matrices Hj defined by is valid. Thus, the matrices to which the 
method described here can be applied are (with some unimportant restric- 
tions) the same as the matrices, to which the Moiseev's method is applicable. 

Let us formulate one important consequence of the property (J7J). 

Proposition 1 If property |7|) is fulfilled then there exists rational matrix 
B(k) commuting with all matrices Hj(k). 

The proof of the proposition is as follows. Represent Hi in the form 

H 1 (k) = P(k)F 1 (k)p-\k), (9) 

where P(k) is the matrix composed of the eigenvectors of H\ normalized, 
say, by making the first component of each vector equal to 1. Respectively, 
F\{k) is a diagonal matrix composed of scalar functions fi(k), . . . , /jv(^)- 

It is known that if two matrices commute then normalized eigenvectors 
of the matrices coincide [Hj. Matrix P is algebraic, so it is single- valued on 
some Riemann surface. Since the values of Hj(k) taken on different sheets 
(with the same k) commute, we can conclude that when a branch point of P 
is bypassed the columns of P are just permuted. 
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Construct matrix B in the form 



B(k) = P(k)D(k)P-\k), 



(10) 



where D(k) is a diagonal matrix with the scalar functions hi(k), . . . , hw(k) on 
the diagonal. Let the functions h m be branches of some algebraic function 
h having the same branch points as P. Moreover, let the values h m be 
permuted the same way as the columns of P when the branch points are 
bypassed. Then the function B(k) is single- valued, and therefore rational. 
A proper choice of the functions h m is as follows: 



where n (k) is an arbitrary set of scalar rational functions (provided none of 
h m is identically zero). Values of P n ,m are elements of matrix P. 

Since B commutes with Hi and all matrices Hj have distinct eigenvalues 
almost everywhere, matrix B commutes with every Hj. Note that all other 
matrices Hj can be represented as 



Due to arbitrariness of the choice of j3 n (k) one can make matrix B having 
simple poles and tending to / as \k\ — > oo. The matrix B possessing all these 
properties plays an important role below. 

4 A family of Riemann— Hilbert problems and 
derivation of ODE1 

4.1 Family of Riemann— Hilbert problems 

We have reduced the matrix factorization problem to finding the function 
U(k) obeying equations (jSJ) on the cuts Tj. To solve this problem we use the 
idea described in detail in [15]. Namely, we are fixing the functions Hj(k), 
defined and continuous on contours Tj (and regular in the strips Q + kj), and 
introduce truncated contours 



N 




(11) 



n=l 



Hj(k) = P{k)F j (k)p- 1 {k). 



(12) 



Tj(b) = (kj + b, kj + ioo), 



7 



where 6 is an imaginary number 6 G (0, zoo). Consider a family of problems 
for the function {7(6, k) set by the relations 

{7(6, k + ) = U(b, k-)Hj(k), k e Tj(b). (13) 

We assume that for each b the matrix function U(b, k) is single-valued, con- 
tinuous, and free of zeros of determinant on the plane of k cut along the 
contours Tj(b). It tends to I as \k\ — » oo. We assume also that the behavior 
of U(b, k) at the points kj + b is derived by continuity from the conditions 
formulated for large Im[6] (see below). Obviously, 

U(k) = {7(0, k). (14) 

The main idea of the method is to study the behavior of U(b, k) as a function 
of b. 

Embedding of U(k) into a family {7(6, k) enables us to define behavior of 
U(k) at k — kj in the most natural way. Expand matrices Hj in the form 
( I12p . The leading term of {7(6, k) near the point k = kj + b can be written 
in the form 

{7(6, k) « K 3 (kj + b)(k- (kj + &))M^+&))/(2™)#-i(fc. + 6 ) ( 15 ) 

where Fj are taken from ([2]), and i£j are some non-singular matrices. The 
branch of the logarithm should be fixed as follows. For 6 — > ioo the matrices 
Fj(kj + b) tend to I. For these values we choose the branch of logarithm close 
to zero matrix. Then, for other values of 6 we choose the branch of logarithm 
by continuity. Such a choice enables us to avoid discussing partial indices of 
the initial Riemann-Hilbert problem. 



4.2 Form of ODE1 for a single cut 

Let the number of cuts p be equal to 1, i.e. let there exists only one cut 
Tj. This corresponds to a function G(k) having a single branch point in the 
positive half-plane. This case has been studied in [15]. Here we formulate 
the main theorem of [15] with a short proof. 

For this, it is necessary to introduce a notation of the ordered exponential 
(the term comes from quantum mechanics). Namely, let 7 be a contour 
(directed one) connecting the points T\ and r 2 {j\ is a starting point), and 
let C(t) be a iV x N matrix defined on 7. Consider a matrix equation 

4-X(t) = C(t)X(t) (16) 
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taken with the initial condition X(j~i) = I. Solve this equation along contour 
7 and define the value X(r 2 ). By definition, 



0E 7 [C(r)dr] = X{t 2 ). 



(17) 



This notation is just a convenient way to refer to a solution of an ordinary 
differential equation. 

Theorem 1 a) There exists N x N matrix Si(fe) analytical in the strip Q, 
such that U (b, k) obeys an ordinary differential equation ( ODE1 ) 



(18) 



The initial condition for this equation is as follows: 

U(ioo, k) = I. (19) 
b) Let there exist a N x N matrix si(r) analytic in Q and such that 



OE, 



k - (r + h) 



-dr 



Hi(k) 



(20) 



fork G (k±, ki+ioo). Contour^ is a concatenation of^ + and r y (see Fig.\^). 
Then solution U(b, k) is given by the formula 



U(b, k) = OE. 



76 



si(r) 



_k-(r + ki 
where contour 7& goes from ioo to b along Fx. 



-dr 



(21) 



lm[b] 



r 



Figure 2: Contours 7 + and 7" 
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Let us outline the proof. Consider part a). Consider the function 



S(b,k) = ^±U-\b,k). (22) 

This function is analytic in the plane cut along (kj + b, kj + zoo). Consider 
its behavior at the cut. Note that the coefficient Hi does not depend on b. 
Thus, 

Using this relation with ([TBI we conclude that S(b,k + ) = S(b,k~), and 
therefore the function is single-valued. According to the condition at infinity 
for U, S should decay as \k\ — > oo. The only singularity of S in the finite 
part of the /c-plane is k = k\ + b. The leading term of the singularity is given 
by fTTS]) . According to this, S has a simple pole at k — k± + b, and 

S(b, k) = - /^(fci + b) hgiF.ih + b)) K^{k x + b). (23) 

Finally, 

si(6) = -^-Ki{h + b) hg^ih + b)) K^(h + b). (24) 

The choice of branch of the logarithm has been discussed in the previous 
subsection. 

Analyticity of the coefficient si in Q + ki follows from the fact that the cut 
on which functional equation (Tl3|) is set can be deformed (without changing 
its starting point) arbitrarily within Q + k\, and the solution remains the 
same while the contour changes. 

The initial condition (1T0]) follows from general properties of the Riemann- 
Hilbert problem [9]. 

Consider part b) of the theorem. Let us show that (12T1) is a solution of 
the problem f fl3|) . Note that due to analyticity of s\ the contour 7 in (TSUI) 
can be deformed provided that it does not cross the singularity r = k — k\ 
of the coefficient. Let be k G ^i{b), and thus lm[k] > Im[A;i + b]. Deform 
contour 7 into 7 + 6: 

0E 4^4r^H = " lW (25) 
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According to general properties of the ordinary differential equations and the 
ordered exponential notations [To] . 



OE 



■y+b 



k-(r + hi) 



-dr 



OE 



6+7" 



k- (r + h 
Note that according to (|2"T|) 

OE 



-dr 



OE, 



k- (r + ki 



-dr 



si(r) 



OE 



6+7" 



fc - (r + fei 



-dr 



k — (r + h 
Thus, (|2SJ) is equivalent to ffT51) . 



■dr 



£7(6, fc H 
£7(6, jfe" 



(26) 



4.3 Form of ODE1 for several cuts 

If there are p > 1 branch points fcj (and, thus, several cuts Fj) Theorem 1 
can be modified, while the reasoning remains basically the same. Coefficient 
S from fl22|) can be proven to be single-valued and decaying, but it should 
have p simple poles k = kj + b. Therefore equation f fl8|) has form 







-k -(kj + b) 



U(b, k). 



(27) 



with p unknown matrices Sj(b) analytical in Q. The initial conditions are the 
same as for one cut (i.e. (j!9p ). The form of the solution follows from (123 
and (USD: 



£/(6, fc) = OE 



ib 



.i=i 



(r + %) 



-dr 



(28) 



A generalization of ( l24j) has form 

1 



Si (6) 



2vri 



+ b) log^ffcj + 6)] Kj 1 (k j + 6). 



(29) 
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for some unknown matrices Kj and known (up to transmutations) diagonal 
matrices Fj. 

Finally, condition (120]) should be rewritten as 

= H m (k), keT m + b, m = l...p. (30) 

5 ODE2 

5.1 Derivation of ODE2 

Formula (12~TT) (or (1281) ) cannot be used immediately to find solution U{k) = 
U(0, k) since matrix functions Sj(b) are unknown. Thus, before finding U one 
should find Sj somehow. In [15] it has been proposed to use equation f l20|) to 
find s\. A numerical procedure has been proposed and tested. Application 
of this procedure does not require branch-commutativeness, so the method 
is potentially applicable to a much wider class of problem than Moiseev's 
class. However, this procedure is rather sophisticated and it does not reveal 
the mathematical nature of the solution. Here we are proposing another 
technique reducing the determination of Sj to solving a (nonlinear) ordinary 
differential equation. Unfortunately, the new technique is applicable only to 
Riemann-Hilbert problems obeying relations ([7J). 

The key idea of the new method is to use matrix B{k) defined by ( TTOT) . 
Namely, we construct a rational matrix B(k) commuting with Hj(k) on all 
their sheets, behaving as B(k) — > I as \k\ — > oo and having only simple 
poles. Obviously, such matrix can be constructed by using the arbitrariness 
of the rational functions (3 n {k). Let the poles of B{k) be located at the points 
k — pi, I — 1 . . . d. 

Consider function 

V(b,k) = U(b,k)B(k). (31) 

Note that 

^P^'(M) = ^W*) S S(t,k). (32) 

Thus, V obeys ODE1 ( 1271) with the same coefficient as U. 

The key property of function V is expressed by the following proposition. 



OE„ 



E 

,i=l 



SAT 



kj) 



dr 
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Proposition 2 There exists function Rib, k), which is rational as a function 
of k for each b, such that 



V(b,k) = R(b,k)U(b,k). (33) 

Construct function R as follows: 

R(b, k) = V{b, k) £7~ 1 (6, k) = U{b, k) B{k) U-\b, k) (34) 

Consider the behavior of V(b, k) on the cuts Tj(b). Since B commutes with 
all Hj, 

R{b, k+) = £7(6, AT) Hj(k) B(k) Hj\k) £7~ 1 (6, k~) = 
£7(6, fc") B(k) £7- x (6, k~) = R(b, k~), k e Tj(k). 

Thus, for each b function R(b, k) is a single-valued function of k. At infinity 
R(b,k) —> I. Obviously, R can only have simple poles at k = p\. Due to 
Liouville's theorem, R(b, k) should be a rational function of k. Moreover, one 
can conclude that R has form 

R(b,k)=I + J2r^-, (35) 

where r^(6) are some N x N matrix functions of b defined in Q. 

Construct the coefficient of ODE1 for V using representation (1331) : 

ab ob ob ob 

Comparing ( 1361) with ( |32|) . conclude that 

dR(b, k) 



db 



SR-RS = [S,R}. (37) 



Equation ( |371) is the global form of ODE2. One can easily see that this 
equation describes the evolution of R but from the first glance it is not clear 
how it can describe the evolution of S. However, we possess some additional 
information about R and S (namely, both functions are rational with respect 
to k). This information is enough to transform fl37|) into a local form, which 
is a closed set of ordinary differential equations describing the evolution of 
R and S. 
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Substitute ( 157)) and 

*ft*) = EizfW (38) 

into (157)) . Expand the right-hand side and left-hand side of (157)) as a sum 
of simple fractions. Taking into account that 

1 1 



k — (kj + b) k — pi (kj + b) — pi \k — (kj + b) k — p t 

and considering the terms with each denominator separately, obtain equa- 
tions 

f ^MM, i = 1 ... d (39) 

and 

= i = l...P (40) 

j-f Pi - (% + 6) 

System (1591) . (|4Q|) does not form a closed system of ordinary differential 
equations for finding the unknown matrices Sj(b), rj(b). To make the system 
closed, consider ( 140)) together with ( 129)) . Formulate the problem of finding 
of matrices Sj provided that matrices r\ are known. Equations (129)) provide 
information about the eigenvalues of Sj, while ( l4"0l) provide information about 
the eigenvectors of Sj. Namely, the eigenvalues of Sj(b) are equal to the 
diagonal elements of 

F 3 (b) = -^-log(F J (k J +b)). 

where Fj(kj + b) is a (known) diagonal matrix composed of the eigenvalues 
of Hj{kj + b) (see (EE2J)). According to (140)) . the eigenvectors of Sj(&) coincide 
with the eigenvectors of the matrix 

R{b,kj + b)-^ 



^{kj+b)-pi 

Define function F(X, Y) producing a matrix, whose eigenvalues coincide 
with the eigenvalues of X, and the eigenvectors coincide with the eigenvec- 
tors of Y (provided all eigenvalues of Y are distinct). The function T is 
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defined ambiguously since the mapping between the eigenvalues of X and 
eigenvectors of Y is not defined. I.e. J 7 is defined up to a permutation of 
order N. If this ambiguity is eliminated in a correct way, 

8 j = ?(F j (b),R(b,k j + b)). (41) 

Equations f l39j) together with f|4Tl) form a closed system of equations for 
finding r\ and Sj. This system is non-linear. The system ( 139]) . (j4"TT) will be 
called the ODE2 (in the local form). Derivation of the ODE2 is the main 
result of this paper. This result can be formulated in the form of the following 
theorem. 



Theorem 2 Let there be a family of Riemann-Hilbert problems |3]) obey- 
ing the restrictions posed above, including the commutativity restrictions |?|). 
The ODE1 for this family has the notation of (2l\ ). Then there exist such 
matrices ri(b) and such a choice of the function J 7 that matrices Sj(b), ri(b) 
obey the system (dPP, (R^P- 



5.2 Initial conditions for ODE2 and choice of func- 
tion T 

To make a numerical solution of ODE2 possible one should define the initial 
conditions and eliminate the ambiguity of defining the function J 7 . Since 
U(b, k) — )■ / as b — )■ zoo, one can conclude that 



where of course 



Thus, if 



R(ioo,k) = B(k), (42) 



R(ioo,k) = lim R(b,k). 

b—tioo 



i=i Hl 



for some matrices U (which are assumed to be known) then 

rj(ioo) = t h (44) 
These relations play the role of initial conditions for the ODE2. 
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To eliminate the ambiguity of definition of function J 7 , we need to estab- 
lish a correspondence between the eigenvectors of the matrix R(b, kj + b) and 
the diagonal element of the (diagonal) matrix Fj(b). Again, consider large 
values of Im[6]. For large imaginary b the values of the coefficients Hj(k) 
approximately commute with the common factor B(b). Therefore, the solu- 
tion U(b, k) near the points k = kj + b approximately commutes with B(b) 
or (which is the sam in asymptotic sense) with B(kj + b). Thus, according 

to (USD, 

R(b,kj + b) « B(kj + b). 

Using this relation, one can establish a natural correspondence between the 
eigenvectors of B(kj + b) and R(b, kj + b). Then, the eigenvectors of Bikj + b) 
are by construction the eigenvectors of Hj{kj + b). Thus, it is possible to 
establish a natural correspondence between the eigenvectors of B(kj + b) and 
Hj(kj + b). Finally, this gives correspondence between the diagonal elements 
of Fj{kj + b) and the eigenvectors of R{b, kj + b). 

Thus, function J 7 can be defined without ambiguity for large Im[fe] and 
for other b it can be defined by continuity. 

5.3 Invariance of ODE2 with respect to the choice of 

B(k) 

The choice of the factor B(k) is not unique. Namely, if B(k) obeys all 
restrictions then a combination 

N-l 

B'{k) = Y,9 m {k)B m {k) (45) 

m=0 

with rational scalar functions g m {k) also can be used as B, provided that 
B'{k) — > I as \k\ — > oo and B' has only simple poles. The form of ODE2 
changes when B is substituted by B' . Let us show that this substitution 
does not change the solution Sj{b). For this we remind that the system (|39|) . 
(14"U|) is equivalent to (1371) . The invariance of Sj is established by the following 
proposition. 

Proposition 3 Let R be defined by (34\ ), and 

R'(b, k) = U{b, k) B'{k) U-\b, k), (46) 
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where B' is defined by (JW . Let [Sty be valid for some matrix S . Then 

dR'(b,k) 



Ob 

First, note that it follows from (146 p that 

N-l 



[S,R']. (47) 



R\k) = Y,9 m {k)R m {k) (48) 



m=0 

Due to formal linearity of (14"T|) . it is sufficient to prove that 

dR m (b, k) 



db 

This can be easily proved by induction. 



[S,R m ]. (49) 



6 Examples 

6.1 Description of the numerical procedure 

The numerical procedure straightforwardly follows from Theorem 1 and The- 
orem 2. Assume that matrices Hj(k) are known explicitly, and let the matrix 
B(k) be constructed and represented in the form f|43l) . 

First, ODE2 is solved along the positive imaginary axis of b from zoo to 0. 
In practice, ODE2 is solved not from b = ioo, but from b = iL, where L is 
a large number playing the role of infinity. At the "infinite" point b = iL 
initial condition for ODE2 are set in the form of 

n{iL) = U. (50) 

At the point b = iL function T is constructed without ambiguity as fol- 
lows. According to the argument above and according to (1501) . for numerical 
solution 

R(iL,kj +iL) = B(kj + iL). (51) 
Represent B(kj + iL) in the form ffTOl) . i.e. 

B(kj + iL) = P*D*p-\ (52) 
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where is a diagonal matrix. Function Hj(kj + iL) can be represented in 
the form (1121) . i.e. 

P~ 1 H j (k j +ib)P, 
should be a diagonal matrix. According to (jUJ, 

Sj (iL) = ~P, log^Hjik, + ib) P*)P-\ (53) 
ziri 

where the branch of logarithm close to zero is taken. This procedure defines 
Sj(iL) in a unique way. 

Then ODE2, i.e. the system ( 139|) . (jHJ) is solved numerically, say by 
Runge-Kutta method, from b = iL to b = 0. On each step function T 
is chosen such that new values Sj(b — i8) are close to old values Sj(b), i.e. 
such that Sj(b) are continuous. As the result of this procedure, the matrices 
Sj are found at points covering the segment (iL, 0) densely enough. 

Next, ODE1 is solved to find U(b, k). A set of points k = z n at which the 
solution U (k) will be found is selected. The initial conditions have form 

U(iL,z n ) = I. (54) 

Equation (128"]) is solved from b = iL to b = for the values U(b, z n ) along the 
imaginary axis of b say by Runge-Kutta method. As the result, the solution 
U(k) = {7(0, k) becomes known at the points k = z n . 

One can see that the numerical procedure is rather simple. It consists 
of two solutions of ordinary differential equations. If the segment (iL, 0) is 
split into Nf, steps, and if there are points in the set z n , then the first 
step takes ~ N b operations, and the second step takes ~ N k N p operations. 
This makes difference with results of [15] where the first step takes ~ N% 
operations. 

6.2 Khrapkov's case 

It is important to show that the proposed technique is equivalent to the 
known method in the simplest commutative case, namely in Khrapkov's case 
[5]. Consider as an example a family of Riemann-Hilbert problems set on 
Ti(b) = (k\ + b, k\ + zoo) with the coefficient 

H 1 (k) = g (k)I + g 1 (k)A(k), (55) 
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where 

A <*> = ( I -i 



go and gi are some algebraic functions such that go(k) — > 1 as |&| — >■ oo, g>i 
tends to zero as \k\ — > oo not slower than l/|fc| 2 . 

A traditional solution of this problem is as follows. First, a solution 
£7(6, /c) is constructed by the formula [5] 

£7(6, fc) = exp(0 [cosh (v^i?) / + sinh (yW)v) ~^^\ ' ( 56 ) 



4>(k) = k 2 + 1, 



fcl+ioo fei+ioo 

^b,k) = - J iM-dr, 77(6, /c) = - | ^dr, (57) 
= log {g 2 (k) - <P(k)9 2 i(k)) , (58) 



Aniy/m \g (k) - gi {k)y/m I 



(59) 



Solution £7 obeys all conditions except the condition £7 — > I at infinity. 
Instead, for a fixed 6 

£7(6, fc) -> cosh(C(6))/ + sinh(C(6)) ^ J J J = Q(6) as ->• 00, (60) 

C(6) = - I i?(r)dT. (61) 

ki+b 

Thus, one has to "correct" the behavior of £7 by a left multiplication: 

U(b,k) = Q- 1 {b)U{b,k). (62) 

Let us consider the same problem from the point of view of the proposed 
method. One can check directly [15] that the auxiliary solution obeys ODE1 
in a slightly modified form: 

(63) 
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Similarly, it can be checked that Q(b) obeys equation 



(1Q i7(6 + *i)f ? n 1 Q(b). (CI) 



db 17 V 1 



Construct ODE1 for [/. According to ( 1621) the coefficient of this equation 
is equal to 

S(b, k) = ^ = Q- 1 (s - ^Q" 1 ) Q = Q- 1 (s-vib + h)^ Ij) Q. 

(65) 

One can see that the coefficient has form of (fTB"l) . i.e. for a fixed b it is a 
rational function of k having a simple pole at k = b + k\ and decaying at 
infinity. 

Now consider ODE2. Select a function B(k) commuting with both branches 
of Hi(k), having only simple poles and tending to I at infinity. For example 
one can choose 

^) = / + F^I AW (66) 

with simple poles at k — ±1. Define R = UBU^ 1 . One can see that B 
commutes with U, and thus R(b, k) = B(k). Note that 

[R, S]=0, = 0. (67) 

Define R as (134j) . It can be expressed as 

R(b, k) = Q-\b)R(b, k)Q(b) = Q-\b)B(k)Q(b) (68) 

Taking into account (1671) and ( 1651) it is easy to show that equation (1371) i.e. 
ODE2 in the global form is valid for Khrapkov's matrix. 

Let us write down ODE2 in the local form (for demonstration purposes). 
Represent A(k) in the form 



A{k) = Vk^+ip(k) n _° x ^ P- 1 



(*), (69) 



20 



Similarly, 

H x {k) = P{k)F x {k)P-\k), 



p m_ f g (k) + VW+i 9l (k) o \ 

llj "V g (k) - y/W+l 9l (k) ) [{L) 

Since R has two poles (pi = 1, p 2 = —1), we need a system of equations 
describing evolution of three matrices: ri(fe), r 2 (&), and si(6). According to 
fl39|) . first two equations have form 



dn(6) [si(6),n(fo)] dr 2 (6) _ [si{b),r 2 {b)\ 



(72) 



rf6 l-^x + fo)' db -1 -(A4 + &)' 

The third equation has form of (HT|) : 

, (6) = _^( w+h)) , j ^ + _^_ ) ). (73) 

Initial conditions for (1721) should be taken in the form (jSJ). For this, 
matrix B should be represented as a sum of simple fractions. As the result, 
we get 

r^) = \{ \ \ ), r 2 (ioo) = \ ( - 1 ^ (74) 

Function F(X, Y) is implemented as follows. Let X be a diagonal matrix. 
Matrix Y is represented in the form Y = Y\ Y 2 Y x ~ l numerically or analytically 
(Y 2 should be a diagonal matrix). The result is formed as 

F(X,Y) = Y 1 XY 1 -\ (75) 



or 

JF{X,Y) = Y 1 X'Y 2 l (76) 

where X' is a matrix, whose diagonal elements are interchanged. The choice 
between these two forms is made by the following rule. For the point b = iL 
where conditions (1741 are set matrix Y has two eigenvectors, one of which is 
close to 



and another one is close to 
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(see fl70p ). These eigenvectors are columns of Y. If the first column cor- 
responds to the vector of the first type, then form (j75p is chosen at this 
point. Otherwise, form (176]) should be chosen. At each new step function T 
is chosen to be approximately continuous. 



6.3 Factorization of Antipov's matrix 

Here we consider a more sophisticated (but also commutative) case previously 
addressed in |3J. Matrix G(k) is as follows: 

G(k)=g Q (k)I + 9l (k)A(k), (77) 

where 

A « - ( %f T ) (78) 

^- m ^$r 3t ' ^mw^rr (79) 

ip(k) = ^k 2 -(l + 0t) 2 (80) 

fi, t, a are some scalar constant physical parameters. Notation (1BU1) means 
that the only branch point in the upper half-plane is ki — 1. 

Matrix (1771) is related to a problem of scattering by a screen composed of a 
rigid half-plane and a flexible perforated sandwich half-plane. The boundary 
conditions for this problem were derived in [16]. The problem was reduced 
to the Wiener-Hopf problem in [4j. 

The problem belongs to the Khrapkov's class. The most important func- 
tion for such problem is 

<j)(k) — A 12 A 2 i - AiiA 22 , (81) 

having the property 

A 2 (k) = (j)(k)I. 

In this case 

(j){k) = k 8 -2fx 4 k 4 + fi s (l + a 2 /r 2 ). (82) 

This function is a polynomial of degree 8. If a direct Khraphov's method 
[5] is applied then a solution grows rapidly (faster than algebraically) at 
infinity. Therefore, the Moiseev's method should be applied. The method has 
been outlined in [4j, however no numerical results have been presented. An 
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application of this method requires finding zeros of Riemann's theta function 
and Weierstrass' kernel quadratures. 

We apply the method developed above to this problem. The following 
values of parameters are taken for computations: 

H = 2, r = 0.25, a = 0.3. 

Apply the Hurd's method. There a single cut Ti in the upper half-plane 
going from k\ — 1 to 1 + zoo. Denote by go(k), gi(k) the values of go(k), 
gi(k) on the right shore of the cut, and by go(k), gi(k) the values on the left 
shore of the cut. Note that these values are different due to the presence of 
the square root if). The coefficient H\(k) describing the multiplicative jump 
on Ti is equal to 

Hl (k) = G(k-)G-\k+) = " *9fti£. ± ~ 9t9 ° )A (83) 

(do) - \9iY<t> 

Then, the function B(k) is chosen. We can take it in the form 

m = I + W) m ' (84) 

where £(k) is a rational function. Since X(k) grows as k 4 , we can take 
as a polynomial of 5th order, namely 

5 

Z(k) = Y[(k- Pl ). (85) 
i=i 

The choice of pi can be done quite arbitrarily. We use the values 

pi = 2 + z, p 2 = 2-i, p 3 = -i, p4 = -l + i, p5 = -l-i. 

The scheme outlined above is implemented. The set of the points of 
interest k = z n belong to the real segment k G (—1,1) (see Fig. [3]). To 
determine them, it is necessary to find the values S\(b) for b G IV These 
values are found by solving ODE2. The result (i.e. the components of the 
matrix £7(0, z n )) is shown in Fig. HI 

Besides finding the values U(0, z n ) we perform a simple control of the 
whole procedure. For this, we find the values U(0, (z' n ) + ) and £7(0, (z' n )~) , 
where the points z' n belong to the cut Ti (see Fig. [3]), the values £7(0, (z' n ) + ) 
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Figure 3: Contours 7 + and 7 



represent the right shore of the cut, and the values U(0, (z' n )~) represent the 
left shore of the cut. To determine the values on the shores we change the 
contour for solving ODE2 slightly. Namely, for the values £7(0, {z' n ) + ) we 
chose contour 7+ in Fig. El and for values U(0, (z' n )~) we chose contour 7". 
After that, we compute the combination t/ _1 (0, (z' n )~)U(0, (z' n ) + )M~ l (z' n ). 
In the ideal case this matrix should be equal to /, therefore its deviation 
from / can be taken as a measure of relative accuracy of the computation. 
It has been found that the relative accuracy of the computation used to be 
of order 10 -4 . 



1.2r 2 




Figure 4: Solution U(k) on the segment k G (—1, 1) 
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7 Conclusion 



A new method for matrix factorization in the commutative (Moiseev's) case 
is developed. The method is numerical, but it is based on two analytical 
properties of the factorization problem. It is applicable to algebraic matrices 
having the property of branch-commutativity, i.e. the matrices, whose values 
corresponding to different sheet over the same affix commute. 

The factorization problem is transformed by Hurd's procedure into a 
Riemann-Hilbert problem on a set of cuts. Then the Riemann-Hilbert prob- 
lem is embedded into a family of Riemann-Hilbert problems indexed by a 
variable b. The solution as a function of b is described by ordinary differential 
equation (ODE1) with an unknown coefficient S. This coefficient is found by 
solving another ordinary differential equation, ODE2. Initial conditions for 
ODE1 and ODE2 are formulated. It is shown that the proposed procedure 
in the Khrapkov's case is equivalent to the standard solution. Moreover, it is 
shown that the new procedure is applicable to Antipov's matrix, and it does 
not lead to Jacobi's inversion problem, which is not easy to implement. 

In more general (non-commutative) cases ODE2 should be replaced by 
an OE-equation described in [15]. 
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